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we extend the so-called Fluctuation Exchange Approximation (FLEX) to study the non-equilibrium 
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FLEX is able to recover all important features of the transition, including the evolution of the 
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systems under non-equilibrium conditions, and, in its rather general form, as formulated here, it 
could find a broad application in molecular electronics calculations. 
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I. INTRODUCTION 

In the past decades, fast and surprising development 
has taken place in the field of molecular electronics. Ex- 
perimentalists succeeded in contacting and gating a va- 
riety of molecules^— , and gained more and more control 
over them. They also managed to fabricate "artificial 
atoms" and molecules from quantum dots, to isolate sin- 
gle electrons on them and manipulate their spin 8 - - — . 

At the same time, theory seems to be legging behind, 
and describing correlated atomic and mesoscopic struc- 
tures under non-equilibrium conditions continues to be a 
challenge for present-day theoretical solid state physics. 
Tremendous effort has been devoted to the development 
of theoretical tools to capture appropriately the transport 
properties and dynamics of these systems jLii~— how- 
ever, with little success. Most methods are uncontrolled 
or work only for rather special models. Under these 
conditions, pcrturbative methods can be of great value: 
Although they are restricted to the regime of weak in- 
teractions, they provide precious theoretical benchmarks 
for more sophisticated though less controlled approxima- 
tions. Furthermore, many experiments are carried out in 
a regime accessible by perturbation theory. 

Theorists typically use the simplest possible models 
such as the (single level) Anderson model or the Kondo 
model to describe correlated behavior in these systems. 
For these simple models it is well-known that pcrturba- 
tive approaches can work rather well in the appropri- 
ate parameter range. In particular, perturbation theory 
in the interaction strength U of the Anderson model is 
known to reproduce the generic structure of the spectral 
functions j22r— although the value of Kondo temperature 
is known to be incorrect Atoms and experimental sys- 
tems are, however, far more complicated than the single 
level Anderson modeL 2 i 28 Typically, magnetic impuri- 
ties contain many electrons on their d or / shells, and 



the orbital structure of these states and the hybridiza- 
tion matrix elements as well as the Hund's rule coupling 
influence quantitatively the corresponding magnetic and 
physical properties. It would thus be important to under- 
stand the limitations of perturbative non-equilibrium ap- 
proaches in multi-orbital systems. Quantum dots, where 
orbital structure can become important under certain 
conditions, offer ideal test grounds in this regard. A 
particular and interesting example is provided by the so- 
called singlet-triplet (ST) transition. 29 ' 30 ! 32 There the oc- 
cupation of two nearby levels (and thereby the spin) of a 
quantum dot with an even number of electrons changes 
due to the presence of Hund's rule coupling. This transi- 
tion has been observed in a number of different systems 
such as vertical 3 - 3 - and lateral quantum dots ; 34 ' 35 carbon 
nanotubesr& or Cqq molecules) 3 - A lot of theoretical ef- 
fort has also been devoted to this transition. In equi- 
librium, the transition can be understood using numer- 
ical renormalization group methods! 31 ' 32 However, our 
understanding of the non-equilibrium situation is rather 
poor: the regime far from the transition could be de- 
scribed through a functional renormalization group (RG) 
approach^ which is, however, not appropriate to de- 
scribe the small bias limit on the triplet side. A slave bo- 
son approach has also been applied relatively successfully 
to describe the somewhat special underscreened case, but 
this approach is rather uncontrolled and is limited to cer- 
tain models^ 

In a previous publication^ we studied the ST transi- 
tion using a simple, perturbative approach, and showed 
that this approach works surprisingly well: It is able to 
capture the physics on both sides of the transition, i.e., 
the two-stage Kondo effect on the triplet side^ . as well as 
the local singlet formation on the singlet side, and the for- 
mation of the corresponding dips in the non-equilibrium 
differential conductance, d//dU. The simple perturba- 
tive approach is, however, not conserving in general ; 39 ' 40 
and furthermore, as mentioned above, it fails to repro- 
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ducc the Kondo temperature^ Therefore, in the present 
work, which should be considered as an extension of our 
previous study, Ref. [H, we go beyond simple perturba- 
tion theory, and study, whether the simplest non-trivial 
conserving approximation, the so-called fluctuation ex- 
change approximation (FLEX) is able to capture the ST 
transition. This method has been extensively applied in 
connection to high-temperature superconductivity ; 4 ^ 42 
and as an impurity solver^ it has also been successfully 
used to combine dynamical mean field theory (DMFT) 
and ab initio techniques It is computationally rel- 
atively cheap, can be extended easily to more than two 
orbitals, and is also able to go beyond perturbation the- 
ory and give a more precise estimate for the Kondo tem- 
perature. 

As we demonstrate, the performance of FLEX is good, 
and it is also able to capture the ST transition. However, 
while it automatically guarantees current conservation, 
its convergence properties seem to be worse than those 
of simple iterated perturbation theory, and it is compu- 
tationally also more demanding. Nevertheless, in spite 
of these weaknesses, FLEX provides a very good option 
to study correlated behavior in nanoscale structures, and 
seems to provide a more accurate estimate for the Kondo 
temperature. 

The paper is organized as follows. In Sees. Ill Al and 
III Bl we introduce the non-equilibrium two-level Ander- 
son model, and describe the fluctuation exchange ap- 
proximation used to solve the non-equilibrium Anderson 
model. In Sec. Ill Cl we show the details of the iteration of 
the full Green's function within the fluctuation exchange 
approximation. In Sec. IIII Al we present the results ob- 
tained for completely symmetrical quantum dots with 
equal level widths, while in Sec. IIII Bl results for dots 
with more generic parameters are discussed. Our con- 
clusions are summarized in Sec. IIVI and some technical 
details are given in the Appendix. 

II. THEORETICAL FRAMEWORK 
A. Model 

Let us start by defining the Hamiltonian we use to 
describe the quantum dot. We divide the Hamiltonian 
into a non- interacting part, Hq, and an interacting part, 
Hint, and write Hq as 

Hq = i?cond + #hyb + #0,dot ■ (1) 

Here the term 

HoAot = E ejd\ a d icr , (2) 

describes the individual levels of an isolated quantum 
dot, and correspondingly, d\ a is the creation operator of 
a dot electron of spin a on level i = ±, with energy 
Si. The other two terms, the conduction electron part, 



-ff con d, and the hybridization, -ffhyb, depend slightly on 
the geometry of the dot. For lateral dots, 

^cond = E^ £,aC£ acr C£ atT , (3) 
ff hyb = E *ai( c U^ + ^ ■ ( 4 ) 

Here £ denotes the energy of a conduction electron mea- 
sured from the (equilibrium) chemical potential of the 
leads, and correspondingly, c^ aa creates a conduction 
electron of spin a in lead a — L,R. In the presence 
of a bias voltage, this energy shifts to £ Q = £ + eV a , with 
V a the electrical potential of lead aJ& The hybridization 
term describes tunneling between the dot level and 
the non-interacting leads, and the parameters t a i char- 
acterize the tunneling amplitude. 

The terms -ff C ond and -f/hyb are slightly different for 
vertical quantum dots or carbon nanotubes. In these 
latter cases, each dot state is associated with a separate 
electron channel in each lead, c^ Q(T — > c^i aa , 

-^cond = £a<4;a<T C C*ao" ! (^) 

H lyb = E t<*i{c\ iaa d ia + h.c.) . (6) 

In this paper, we assume that the occupation of the two 
levels involved in the transition is around a d ia di a ) ~ 
2. Therefore, we write the interaction in an electron- hole 
symmetrical forrn^ 




with U and J denoting the Hubbard interaction and 
the Hund's rule coupling, respectively, and S = 
\ Si a a' diadfjo'dic being the spin of the dot. To carry 
out a systematic perturbation theory, we split the inter- 
action above into a normal ordered term and a level shift, 




We then incorporate the second term in Hq, 




while we treat the normal ordered part 

■■H^:= E \^:fd] a ,dl,,d n ,d ia (11) 

i,j,m,n, 
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as a perturbation. Here the bare interaction vertices, 
n ™ cr can be expressed in terms of U and J, with 
the explicit expressions derived in Ref. H|| The above 
procedure must be contrasted to the one we followed in 
Ref. HH, where the second term of Eq. ([5J has been treated 
through the application of a counterterm procedure. This 
counterterm procedure becomes unnecessary in FLEX, 
which is formulated in terms of the full (dressed) Green's 
functions. 



B. Out of equilibrium fluctuation exchange 
approximation 

To describe the spectral and transport properties of 
the dot, we use a Green's function method. We thereby 
consider the Keldysh Green's functions of the dot elec- 
trons, 

G%:\t t') = i <T^w(*)4 K (*')> , (12) 

with (...) denoting the average with respect to the sta- 
tionary density matrix, T% the time ordering along the 
Keldysh contour, and k and k' = 1,2 the Keldysh indices, 
labeling the upper and lower Keldysh contours. Through- 
out this paper we shall consider the simplest case, where 
the Hamiltonian is spin rotation invariant. In this case, 
the Green's function is spin diagonal, 



G 



(t-t')=8l G\\, (t — t ) 



(13) 



The non-interacting Green's functions, g\ are asso- 
ciated with Hq, and can be determined analytically (see 
Appendix [A] for their explicit form) . They are related to 
the full Green's functions through the Dyson equation, 



G- 1 H = g - 1 H-s( w ) 



(14) 



where we used a matrix notation, m\ — > m, and intro- 
duced the Keldysh self-energy, £. 

Just as in Ref. HH, the knowledge of G enables us to 
compute the current through the dot by using the Meir- 
Wingrccn formula, 



I = — 



h3 . 



die 



(r£-rS)(G<)2(«) + 



(15) 



+ {h{u)v% M w )r*)((G>)?'M - (G<)|H) 

with the lesser and greater Green's functions dchned in 
the usual way in terms of the Keldysh Green's function 
in Eq. JI3J: 



(G > )i 
(G<)i 



.1 -2 



G: 



(16) 
(17) 



The functions / Q (w) = f{u> — eV a ) in Eq. (JT5J) denote the 
shifted Fermi functions in lead a, and the matrices T" 
describe the decay of the dot levels. They arc dchned as 
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2?T N oc tait'a.j i 



(18) 



for lateral quantum dots, while they read as 
(r^-)vcrt = Sij 2n N m \t ai \ 2 , 



(19) 



for vertical dots, with N a and N a i standing for the den- 
sity of states in the leads. We remark that the fac- 
tor N a can be eliminated by incorporating it in the 



7-1/2 



-> tr 



and the fields 



tunneling parameters, t a iN a 

1/2 

Our primary purpose is to determine S (and thus G), 
and use that to compute the non-equilibrium differen- 
tial conductance through the dot. We shall use the so- 
called fluctuation exchange approximation (FLEX) for 
this purpose. FLEX is constructed in terms of a gen- 
erating functional, $ = $[G], defined as a functional 
of the full many-body Green's function, G4£ The self- 
energy and the particle-hole irreducible vertex functions 
are obtained from $ through functional differentiation. 
Although $ is usually not known, one can approximate 
it by a subset of diagrams, and then obtain approxima- 
tions for the self-energy and the vertex functions. As 
shown by Kadanoff and Baym ; 39 ' 40 this construction is 
conserving, i.e., it guarantees that conservation laws are 
respected. Although this approach is mostly used in 
imaginary time, one can quite naturally generalize it to 
the non-equilibrium case discussed here, by simply re- 
placing the imaginary time Green's function in <I> by the 
Keldysh Green's functions. 

In this language, Hartree-Fock theory is just the sim- 
plest conserving approximation, while the next level of 
approximation is provided by FLEX, corresponding to 
the summation of an infinite series of ladder diagrams 
(see Fig. [T|). In Fig. [1] we introduced the Keldysh 
particle-hole vertex, 



£l,<7l,«l /2,CT2,K2 
il.CTl 1-2,0-2 



-S( K l) < 5kiK 2 K 3 K4 f ^ 
p ^35^3 ^2 ,G2 



,&3 li 

,eri h 



(74 
(T 2 



(20) 



with s(k) keeping track of the sign change of the inter- 
action on the Keldysh contour: s(l) = +1 for the upper 



FLEX 





FIG. 1: The $ functional generating the FLEX diagrams. 
The first diagram just generates the Hartree-Fock approxi- 
mation. Heavy lines denote full Green's functions. Squares 
denote the particle- hole vertex, defined in Eq. (f!?U)) . 



4 




FIG. 2: Structure of the particle-hole vertex, F al . Here a, 
P, 7 and 8 denote composite indices as introduced in Eq. (|21[) 



and s(2) = —1 for the lower contour. The structure of 
this particle- hole vertex, T, is shown in Fig. [2] for the 
particular case of Hund's rule coupling and Hubbard in- 
teractions. 

Differentiating the functional $ of Fig. [TJ one obtains 
the self-energy diagrams shown in Fig. [3Ja. We then ob- 
serve that all higher order diagrams contain the ladder 
series, shown in Fig. [31b. Let us therefore introduce the 
composite label, 

(k,(Xi,Ki) -> on , (21) 
and define the particle-hole propagator, n' ) as 

nW^t-Osi'GSCt-OGS^-*)- (22) 

Then the full particle-hole propagator, n, defined by the 
ladder series in Fig. |3Jb. satisfies the following Dyson 
equation: 



(23) 



OO 



Q2.P2 



The integral being just a convolution, this equation can 
be solved in Fourier space. Defining then 



ladder 

-£ 



(24) 

t'), 



=2.^2 



we can sum up all n > 3 order self-energy diagrams. The 
self-energy £ ladder also contains the second order self- 
energy contribution, but with double weight. Therefore, 
the total self-energy can be written as 



J ladder 



E(1)_E(2) 



(25) 



with the first and second order diagrams, E' 1 ) and £^ 2 ) 
defined as 

OO 

= w ;=<w- / ^e^V), (26) 



2tt 

oo a/3 



S (2) f(t-t') (27) 

= -^EE^n^CS:^-*')^^^-*'). 

q./3 ai>ft 

° 2 ,/3 2 

Solving the equations above turns out to be numeri- 
cally rather demanding for two reasons: First, to get a 
good enough time resolution, we have to keep a large 
number of time (frequency) points in the calculations. 
Second, the propagator n has too many indices. In fact, 
even in our simple case, n has 8 4 components. This num- 
ber can be, however, substantially reduced if we exploit 
the SU(2) spin symmetry of the problem. Using simple 
group-theoretical arguments, we can show that the ver- 
tex r assumes a simple form in spin space, and can be 
expressed in terms of a singlet and a triplet component, 



FLEX 



n = 




+ 



=(7l<72 



2 V=s ' =t 










ME, 



E t 
o o 



r v E s - r t 
o 
o 



(28) 



with the four indices ordered as {tt)t4-i4-t>44-}> an d the 



matrices T defined as 

=s,t 



r = 



p tt _ p tt 
= tt = u 

p Tt p tt 
= tt = U 



(29) 
(30) 



+ 



(b) 



+ 



FIG. 3: (a) Series of self-energy diagrams generated from 
$flex- (b) Diagrammatic definition of the full particle-hole 
propagator, n. 



Here each entry is a matrix in the remaining orbital (I) 
and Kcldysh («) labels: (f ^t ts - By the 

same symmetry argument, we can show that the propa- 
gators n(°) and n take on a similar form. Furthermore, 
it is easy to sec that this structure is maintained un- 
der multiplication, where the lower indices of a tensor 
are contracted with the upper indices of another tensor. 
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Therefore the singlet and the triplet components of II can 
be summed up independently: 



n >) = S i°>H 



n (0) M , (3i) 



= s,t = s.t 



with the unit matrix 1 defined as (l)! 3 " 3 ''' 4 " 4 



We can then simply express the spin- 



independent part of Eiadder m terms of II g £ as 



(w);(t) = -S|(£«B,(<)t«)™G|w 



(£.a.wt,)>sw 

P,9 



(32) 



with p and g denoting composite labels, only including 
the orbital and the Keldysh indices, (1,k) ~ > p, q. 



Although for small interaction parameters the conver- 
gence was rather stable, FLEX showed instabilities for 
high interaction parameters, similar to IPT—. These in- 
stabilities could be partially cured by a gradual increase 
of the interaction parameters. With this trick, the range 
of applicability was found to be roughly the same as the 
one found with IPTM- 



III. RESULTS AND DISCUSSION 

Let us now turn to the presentation of the numerical 
results. For simplicity, excepting Subsection 1111 Cl in this 
section we shall focus on a completely symmetrical dot 
with an even (i = +) and an odd (i = — ) level. In this 
case, the tunneling matrix elements satisfy 



±t R ± 



(35) 



C. Details of the FLEX iteration 

The previously defined equations provide a self- 
consistent set of equations, which we then solve itera- 
tively. In O'th order, we approximate the full Green's 
function G by g, 



o . 



(33) 
(34) 



We then start iteration n > 1, by first computing 
G^-^(t) from the Green's function G^^H of the 
previous iteration, by performing a Fast Fourier Trans- 
formation (FFT). Next, we construct {XL { °\) ln ~ 1] (t), ob- 
tain from that QT°')[™ _1 ] (w), and then we solve the 
Dyson equation, Eq. (J3TJ) to get U [ ™~ 1] {u>). From 
that we obtain n tn 7 1] (*) b Y FFT - We can then use 
S [ ™; 1] (0: (II ( s 0) t ) [n " 1] W, and G^-^it) to compute 
^ladder' (£ (1) ) W , and (XjC 2 ))!™], and finally the total self- 
energy, SW(t), through equations Eq. ((Ml), (EH), (EZJ) 
and (f25j) . Finally, we obtain our next estimate, G^(ix>), 
by first computing the Fourier transform, £[™](<j), and 
inverting the Dyson equation, Eq. (fT4")) . This iteration 
procedure is repeated until convergence is reached. 

In the numerical calculations we represented the 
Greens functions using a finite uniform mesh of N fre- 
quency points in the range — 0/2 < uj < 0/2. As men- 
tioned above, the numerics was highly demanding; we 
had to use 2 16 - 2 18 frequency points and w 1000C/ 
to reach convergence. The memory demand of the cal- 
culation was also much higher than that of the iterative 
perturbation theory (IPT) procedure of Ref. [H. With 
the symmetry-based representation of T and II propa- 
gators, however, we managed to reduce the size of them 
substantially and were able to run the calculation on sim- 
ple PC's. 



and the tunnelings can be characterized simply by the 
widths of the levels, 



(36) 



a=L,R 



both for lateral and for vertical dots. Similar to Ref. |38|, 
here we shall focus onto the vicinity of the electron-hole 




0.8- 





A/U=0.0 




A/U=0.3 




AAJ=0.6 




A/U=0.9 




A/U=1.2 






FIG. 4: Total spectral function, pr{oj) = (p+(u>) + p_(w))/2 
for J/U = and Y±/U = 0.785 for different values of level 
splitting, A/U. On the top we show the FLEX results, while 
the bottom shows the IPT results for the same parameters. 
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symmetrical point, e+ = e_ = 0, and assume that the 
two levels are symmetrically positioned, 

e± = ±A/2 . (37) 



A. The case T+ = T_ 

1. Equilibrium, spectral functions 

In this case, for A = J = 0, the three singlet and 
the triplet states of an isolated doubly occupied dot are 
completely degenerate, and an unusual Kondo state is 
formed. 38 ' 47 Turning on A, one separates the singlet state 
with both electrons on state i = — from the rest of the 
states, and destroys the Kondo effect once A becomes 
larger than the Kondo temperature, T£, defined as the 
halfwidth of the central peak for A = J = 0. This tran- 
sition can be observed in the total equilibrium spectral 
functions, 

»=± 

where the retarded Green's function is defined as, 

i,j = ^jl ~ { - r j2 — ~ U i,j ! W y J 



2. Comparison with numerical renormalization group 

Before entering the discussion of the non-equilibrium 
results, it is worth comparing FLEX with other methods 
such as iterative perturbation theory (IPT), or numer- 
ical renormalization group calculations (NRG) , 48 ' 49 the 
latter procedure giving us a benchmark for the equilib- 
rium calculations. Fig. [6] compares the results of these 
three methods for parameters A = 0, T±/U = 0.785, and 
J/U = 0.15. For the NRG calculations we used the open 
access Budapest NRG code<££ To reduce computational 
effort and achieve sufficient accuracy, we made use of the 
spin SU(2) symmetry of the Hamiltonian, as well as the 
U(l) symmetries corresponding to the conservation of 
the total fcrmion numbers in channels i = ±. The com- 
putations were performed with a discretization parame- 
ter, A = 2, and 2400 kept multiplets. The calibration 
of the NRG parameters requires special care, since the 
NRG discretization and iteration procedure renormalizes 
somewhat the bare parameters of the Hamiltonian . 48 ' 51 
We calibrated the level widths T± from the height of 
the numerically calculated spectral functions. The re- 
sults obtained this way were i n g ood agreement with the 
analytical expressions of Ref. [fLU. 

As shown in Fig. [6] the width of the Kondo resonance 
is perfectly captured by FLEX for the above parameters, 
while IPT slightly overestimates the size of the Kondo 
resonance. (As a comparison, in Fig. [5] we also plotted 



with Gf , 



Si 



the time-ordered Green's function. 
In Fig. |H we display pr(ui) for J — for various split- 
tings of the two levels, A, as computed by FLEX and by 
the iterated perturbation theory (IPT) of Ref. [H. The 
splitting of the Kondo resonance is remarkably similar 
in the two figures, however, there are important differ- 
ences, too. First of all, FLEX gives a smaller Kondo 
temperature, and provides a more realistic shape for the 
Kondo resonance both in the absence and in the pres- 
ence of splitting. However, while the Hubbard peaks at 
uj = ±JJ are still visible within the simple perturbative 
calculation, FLEX is unable to capture them correctly. 

Similar conclusions are reached for J ^ with the ex- 
ception that now the splitting of the Kondo resonance is 
shifted to higher values of A (see Fig. [5]). However, in 
this case the central peak has a slightly different inter- 
pretation than for J = 0, since for J > the isolated 
dot would be in a triplet state. As a result, the central 
Kondo resonance at A = can be interpreted as a re- 
sult of a triplet Kondo effect, where the spin S = 1 of 
the dot is screened by the even and the odd conduction 
electron channels. In this triplet state the ground state 
degeneracy of the isolated dot is reduced, and quantum 
fluctuations are therefore somewhat suppressed. As a 
consequence, the Kondo temperature Tk is also reduced, 
and the central peak becomes slightly narrower, but also 
more stable against A ^ 0; in this J > case the split- 
ting of the triplet Kondo resonance occurs roughly when 
A - 2J + T K . 





FIG. 5: Total equilibrium spectral function, pt(^) for J/U = 
0.15 and Y±/U = 0.785 for different values of level splitting, 
A/17, as computed by FLEX (top) and by IPT (bottom). 
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o 

co / (nr/2) 



FIG. 6: Comparison of the spectral functions for A = 0, 
T±/U = 0.785, and J/U = 0.15, as computed by FLEX, IPT, 
and NRG. Clearly, FLEX seems to capture rather accurately 
the width of the central Kondo resonance. 



the shape of the resonance for U = 0.) However, while 
FLEX seems to give a better estimate for the Kondo 
temperature than IPT, IPT seems to capture the high- 
energy features (Hubbard peaks) better - a well-known 
shortcoming of FLEX42 




FIG. 7: Level-resolved (top) and total (bottom) equilib- 
rium spectral functions for J/U = 0.15, T+/U = 1.1 and 
T-/U = 0.785 for different level splittings, A/U, as computed 
by FLEX. The inset in the upper panel shows the normalized 
spectral functions p±(u>) / p±(0) , demonstrating the presence 
of the two different Kondo scales. 



B. The asymmetric case, 7^ T_ 

Let us now turn to the more generic situation, T + 7^ 
r_, and J > 0. In this case, for A « 0, the triplet spin 
on the dot is screened by a two-stage Kondo effect^ 
and the central resonances in the level-projected spectral 
functions, p±(uj), become different due to the presence 
of two different Kondo scales, T^, corresponding to the 
screening in the even and in the odd channels, respec- 
tively. 

In Fig. [7] we show the level-projected as well as the full 
spectral functions, pt(uj), as computed by FLEX for a 
dot with J/U = 0.15, T+/U = 1.1 and T-'/U = 0.785 for 



0.1 
0.075 
0.05 



lateral dl/dV 

"T 



vertical dl/dV 




-2 2 
eV/U 



FIG. 8: Differential conductance, G(V) = dl/dV, in units 
of e 2 /h for lateral (l.h.s.) and vertical (r.h.s.) dots, with 
J/U = 0.15, T+/U = 0.785, and T-/U = 1.1 for different 
level splittings, A/U, as obtained by FLEX. On the triplet 
side (A = 0), the second Kondo scale emerges as a narrow 
dip/sharp resonance in G(V) in the lateral/ vertical arrange- 
ment. For large A's the singlet-triplet splitting gives rise to a 
wide central dip in G(V). The curves reproduce very nicely 
all experimentally observed features, however, the cross-over 
regime is only qualitatively captured in the lateral case. 
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0.5 1.0 1.5 

A/U 



FIG. 9: Linear conductance, for lateral and vertical dots, 
with J/U = 0.15, F+/U = 0.785, and T-/U = 1.1 for differ- 
ent level splittings, A/U, as obtained by FLEX. 

different level splittings, A/U. Unlike for T + = T_, for 
A = the projected spectral functions of the two levels 
are different, p+(co) ^ Nevertheless, they are 

all symmetrical as a consequence of a discrete particle- 
hole symmetry (see Ref. HH). However, this symmetry is 
violated for any A^O, where electron- hole symmetry is 
destroyed even for the total spectral function, pt(oj). The 
difference in the Kondo temperatures is clearly visible in 
the normalized level-projected spectral functions, shown 
in the inset of Fig. [7] 

Similar to the symmetrical case, the Kondo resonances 
are gradually split by a finite A. The splitting of the 
resonances appears even more strikingly in the differ- 
ential conductance, G(V) = dl/dV, as computed from 
Eq. (TTSl) , and shown in Fig. [5] These differential conduc- 
tance curves were obtained by computing G and I(V) 
for each bias voltage V separately, and then carrying out 
a numerical differentiation. 

For a lateral dot at A = 0, i.e. in the two-stage 
Kondo effect regime, the dl/dV curve shows very nicely 
the build-up of the first Kondo resonance ; 34 ' 35 and then 
the appearance of a dip at V = bias. This dip is a 
result of the destructive interference between the two 
Kondo effects, and it appears once the bias voltage be- 
comes so small that it cannot destroy even the narrower 
Kondo resonance of the spectral function. As shown in 
Fig.O increasing A, the linear conductance (i.e., the zero 
bias differential conductance) exhibits a maximum in the 
cross-over regime, in agreement with the experiments. 
However, the bias-dependence of the differential conduc- 
tance in the cross-over regime (dashed lines in Fig. [5]) of 
maximal conductance is not very reliable, and the G(V) 
only shows the general trends observed experimentally, 
i.e., the disappearance of the central dip, and the ap- 
pearance of a state with a single Kondo resonance and a 
perfect G = 2e 2 /h linear conductance. For even larger 
A's, however, the dl/dV curves show very nicely the lin- 
ear splitting of the Kondo resonance. 

In contrast to the lateral case, in a vertical geome- 
try, the second Kondo effect manifests itself as an ad- 
ditional contribution to the conductance, and thus as a 
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FIG. 10: Differential conductance (in units of e 2 /h) of a 
vertical and a lateral dot with J/U = 0.15, T+/U — 1.1, 
Y-/U = 0.785, and A = 0, as obtained by IPT. The curves 
compare quite well with the ones in Fig. [8] 



narrow peak at zero bias for A = 0. In this vertical case, 
the differential conductance curves reproduce the experi- 
mentally observed features even in the cross-over regime: 
the linear conductance is suppressed with increasing A 
(see Fig. [S])^ and the central resonance gets gradually 
broader, until it splits into two side-peaks, corresponding 
to the singlet-triplet excitation energy. 

Finally, for a comparison, in Fig. [TU] we show the 
dl/dV curves at A = 0, as obtained by IPT, for the 
same parameters as the ones used to produce Fig. [8] The 
IPT curves are strikingly similar in structure to the ones 
obtained by FLEX. The most important difference is in 
the width of the central dip/resonance structure, which 
is somewhat narrower in the FLEX calculation, and is 
closer to the real value. 



C. The fully asymmetrical case 

So far, we focused on the case of a completely sym- 
metrical quantum dot, and correspondingly, we assumed 
that one of the states is even while the other state is odd. 
In general, however, quantum dots are not entirely sym- 
metrical. Such asymmetry leads to the suppression of 
the maximal conductance, and for lateral quantum dots 
it may also lead to interference effects^ It is out of the 
scope of the present paper to study such interference ef- 
fects in detail, however, to demonstrate how FLEX works 
in this more general case, let us present here some results. 

In this general case, we can parametrize the tunneling 
to the leads using the angles <f>± G [— tt/2,tt/2] as 

(t±,L,t±, R ) = t± (cos(</5±), sin(<£±)) . (40) 

For an even level, <j> = 7r/4, while for an odd level, cf> = 
-tt/4. 
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FIG. 11: Equilibrium dimensionless spectral functions for 
J/U = 0.15 F+/U = 1.1, T-/U = 0.785, A/U = for 
a symmetrical dot with <f>± = ±7r/4 and for an asymmet- 
rical dot with <j>± = ±7r/3, as computed by FLEX. The inset 
shows the off-diagonal component of the spectral function for 
<j> ± = ±71-/3. 

In Fig. QT] we present the equilibrium spectral func- 
tions, 

^■M^(g«h-g£m), 

for the same level width, T+/U = 1.1, and T_/U = 0.785 
as before, but for a lateral dot with (j>± = ±tt/3. In this 
case left-right symmetry is absent, and pij has offdiag- 
onal components, too. Interference between the states 

± appears as a resonant structure in However, in 

contrast to the components p++ and p , within numer- 
ical accuracy p_| and p |_ integrate to zero according to 

the corresponding spectral sum rule. For A = the dot 
is still electron-hole symmetrical, and the heights of the 
spectral functions at uj = are simply given by 



^•(0) = -(T~ 1 ) ij , 

7T 



(41) 



with Tij — J2 a =L r ^fj t nc Ml relaxation rates (see 



Eq. ((36)) ). as can be checked by an explicit calculation. 

Fig. [T2] shows and compares the differential conduc- 
tance computed for asymmetric vertical and lateral dots 
in the triplet regime (A = 0). The curves are very similar 
to the ones obtained for symmetrical dots, excepting two 
important differences: (a) The conductance of a vertical 
dot does not reach the unitary conductance but goes only 
up to the value 2e 2 /h (sin 2 (20 + ) + sin 2 (20_)) = 3e 2 //i, 
and similarly, the overall conductance of a lateral dot is 
also suppressed, (b) The width of the narrower resonance 
is reduced for a lateral dot. This is due to the fact that 
the smaller eigenvalues of the T matrix are reduced by 
the interference as 



r + r_ cos 2 




and accordingly, the dip corresponding to the narrow 
Kondo resonance becomes also narrower. In contrast, 



FIG. 12: Top: Differential conductance, G(V) = dl/dV, for 
lateral (top) and vertical (bottom) dots with J/U = 0.15, 
T+/U = 0.785, r_/[/ = 1.1 and A/U = 0, as obtained by 
FLEX, for the symmetrical case, <f>± — ±7r/4, and for the 
asymmetrical case, cj>± — ±ir/3. 



the structure of the dl/dV curve remains essentially un- 
altered for a vertical dot, where only the amplitude of 
the signal is reduced. 



IV. CONCLUSIONS 

In the present paper, we developed a general non- 
equilibrium fluctuation exchange approximation (FLEX) 
formalism. We tested the performance of this approach 
on the singlet-triplet transition of a dot with two single- 
particle levels, driven by a competition between the 
Hund's rule coupling and the Kondo screening. This 
transition exhibits several correlation-induced features, 
which are typically rather difficult to capture. On the 
triplet side of the transition a Kondo state develops with 
two different Kondo scales, while on the other side of the 
transition the triplet excitation appears as a pseudogap 
feature. Finally, in the cross-over region an exotic Kondo 
state appears, and for a lateral dot the linear conductance 
shows a broad resonance. 

Remarkably, within its range of convergence, FLEX 
was able to capture all these features, excepting the 
Hubbard peaks, which are rather poorly represented by 
FLEX. Nevertheless, the low energy features and the 
dl/dV curves show behaviors remarkably close to the 
experimentally observed ones. In our earlier studies, we 
applied simple (iterative) perturbation theory (IPT) to 
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describe the singlet-triplet transition. FLEX has some 
clear advantages, but also disadvantages with respect to 
IPT. On the one hand, it produces apparently more re- 
alistic curves in the small bias region than IPT, and - as 
our comparison with NRG calculations confirms - it cap- 
tures the Kondo temperature as well as the Kondo effect- 
related structures better there. In addition, it is a gener- 
ically conserving approximation, and it scales rather well 
with the number of orbitals. All these properties make 
FLEX a viable route to incorporate strong correlation ef- 
fects in molecular electronics calculations. On the other 
hand, FLEX is computationally much more demanding. 
In fact, in this work we had to exploit symmetries to 
reduce the computational effort. This is, of course, not 
a major obstacle if one has access to supercomputers or 
efficient computer clusters, and we believe that the nu- 
merical efficiency can most likely be further improved. 

Finally, let us comment on the version of FLEX we 
used here. In the present paper, we used a generating 
^-functional, which only incorporates electron-hole bub- 
ble series. FLEX can, however, be extended to include 
fluctuations in the Cooper channel, too. This may be 
important in cases, where attractive interactions appear 
in some scattering channels. In particular, such an ex- 
tension of FLEX may be necessary to describe transport 
through superconducting grains. The generalization is 
relatively straightforward, however, it is certainly beyond 
the scope of the present work, which solely focused on the 
demonstration of FLEX as an efficient non-equilibrium 
impurity solver. 
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Appendix A: The hybridized Green's function, g 

For completeness, let us give here the elements of 
g _1 (w). Restricting ourselves to the spin symmetrical 

case, g^Z* = 9^1k ■ The elements of cT 1 ^ dif- 
fer for lateral and vertical dots. For lateral dots, they are 
given by 

(Slat)-;' = 

- Y, N a t* ai t aj Af(u;) , (Al) 

aeL.R 

with s(k) the Keldysh sign defined in the main text, and 
hybridization parameters A£ K (u) defined as 

AiV) = 7ri(2/ a (w) - 1) , (A2) 

Ai 2 (uj) = -2mf a (u), (A3) 

A^V) = -27rf(/ a (w) - 1) , (A4) 

Af(os) = 7ri(2/ a (a;) - 1) , (A5) 

with f a (ui) = f(uj — eV a ) the shifted Fermi function. For 

vertical dots, on the other hand, g~ x Z is diagonal in i 
and j, 

(<?vert);«' = * Ufa) (w - ii) <£' 

- S{ J2 A^M 2 Ar'M- (A6) 

a£L,R 
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